# coding: utf8
# Graham scan 算法，根据2维平面上n个点，求包括这n个点的凸形区域

import math


# points中纵坐标最小点的索引，若有多个则返回横坐标最小的点
def get_bottom_point(points):
    min_index = 0
    n = len(points)
    for i in range(0, n):
        if points[i][1] < points[min_index][1] or (
                points[i][1] == points[min_index][1] and points[i][0] < points[min_index][0]):
            min_index = i
    return min_index


# 按照与中心点的极角进行排序，余弦，center_point: 中心点
def sort_polar_angle_cos(points, center_point):
    n = len(points)
    cos_value = []
    rank = []
    norm_list = []
    for i in range(0, n):
        point_ = points[i]
        point = [point_[0] - center_point[0], point_[1] - center_point[1]]
        rank.append(i)
        norm_value = math.sqrt(point[0] * point[0] + point[1] * point[1])
        norm_list.append(norm_value)
        if norm_value == 0:
            cos_value.append(1)
        else:
            cos_value.append(point[0] / norm_value)

    for i in range(0, n - 1):
        index = i + 1
        while index > 0:
            if cos_value[index] > cos_value[index - 1] or (
                    cos_value[index] == cos_value[index - 1]
                    and norm_list[index] > norm_list[index - 1]):
                temp = cos_value[index]
                temp_rank = rank[index]
                temp_norm = norm_list[index]
                cos_value[index] = cos_value[index - 1]
                rank[index] = rank[index - 1]
                norm_list[index] = norm_list[index - 1]
                cos_value[index - 1] = temp
                rank[index - 1] = temp_rank
                norm_list[index - 1] = temp_norm
                index = index - 1
            else:
                break
    sorted_points = []
    for i in rank:
        sorted_points.append(points[i])

    return sorted_points


# 返回向量与向量[1, 0]之间的夹角-从[1, 0]沿逆时针方向旋转多少度能到达这个向量
def vector_angle(vector):
    norm_ = math.sqrt(vector[0] * vector[0] + vector[1] * vector[1])
    if norm_ == 0:
        return 0

    angle = math.acos(vector[0] / norm_)
    if vector[1] >= 0:
        return angle
    else:
        return 2 * math.pi - angle


# 两向量的叉乘
def coss_multi(v1, v2):
    return v1[0] * v2[1] - v1[1] * v2[0]


# 求多个位置点形成的凸形区域
# Input: points 凸形区域所有点，顶点形式为 [lng, lat]
# Output: 凸形区域边界点
def graham_scan(points):
    bottom_index = get_bottom_point(points)
    bottom_point = points.pop(bottom_index)
    sorted_points = sort_polar_angle_cos(points, bottom_point)

    m = len(sorted_points)

    if m < 2:
        # print("点的数量过少，无法构成凸包")
        if m == 0:
            return [bottom_point]
        if m == 1:
            return [bottom_point, sorted_points[0]]

    stack = []
    stack.append(bottom_point)
    stack.append(sorted_points[0])
    stack.append(sorted_points[1])
    # print('当前stack', stack)

    for i in range(2, m):
        length = len(stack)
        top = stack[length - 1]
        next_top = stack[length - 2]
        v1 = [sorted_points[i][0] - next_top[0], sorted_points[i][1] - next_top[1]]
        v2 = [top[0] - next_top[0], top[1] - next_top[1]]

        while coss_multi(v1, v2) >= 0:
            if length < 3:     # 加上这两行代码之后，数据量很大时不会再报错
                break          # 加上这两行代码之后，数据量很大时不会再报错
            stack.pop()
            length = len(stack)
            top = stack[length - 1]
            next_top = stack[length - 2]
            v1 = [sorted_points[i][0] - next_top[0], sorted_points[i][1] - next_top[1]]
            v2 = [top[0] - next_top[0], top[1] - next_top[1]]
        stack.append(sorted_points[i])

    return stack


# 求凸形区域中心点
# Input: points 凸形区域所有顶点，顶点形式为 [lng, lat]
# Output: 凸形区域中心点
def convex_center(points):
    # 求平均数  同时角度弧度转化 得到中心点
    x = 0       # lon
    y = 0       # lat
    z = 0
    lenth = len(points)
    for lon, lat in points:
        lon = math.radians(float(lon))
        lat = math.radians(float(lat))
        x += math.cos(lat) * math.cos(lon)
        y += math.cos(lat) * math.sin(lon)
        z += math.sin(lat)
        x = float(x / lenth)
        y = float(y / lenth)
        z = float(z / lenth)
    center_coors = [math.degrees(math.atan2(y, x)), math.degrees(math.atan2(z, math.sqrt(x * x + y * y)))]
    return center_coors


# 根据两个点的经纬度，计算两点间的距离
# Input：点1纬度 lat1，点1经度 lng1，点2纬度 lat2，点2经度 lng2
# Output：位置点间的距离 km
def coors_distance(lat1, lng1, lat2, lng2):
    earth_radius = 6378.137
    rad_lat1 = lat1 * math.pi / 180.0
    rad_lat2 = lat2 * math.pi / 180.0
    a = rad_lat1 - rad_lat2
    b = (lng1 * math.pi / 180.0) - (lng2 * math.pi / 180.0)
    s = 2 * math.asin(math.sqrt(math.pow(math.sin(a / 2), 2) + math.cos(rad_lat1) * math.cos(rad_lat2) * math.pow(math.sin(b / 2), 2)))
    s = s * earth_radius
    return s           # km


# 利用海伦公式计算凸多边形的面积
# Input: points 凸形区域所有顶点，顶点形式为 [lng, lat]
# Output: 凸多边形面积 km2
def heron_area(points):
    # 基于海伦公式计算多边形面积
    area = 0
    if len(points) < 3:
        return area

    if points[0] == points[len(points) - 1]:
        pb = ((points[len(points) - 2][0]+points[0][0])/2, (points[len(points) - 2][1]+points[0][1])/2)   # 基准点选为第一个点和最后一个点连线边上的中点
    else:
        pb = ((points[len(points) - 1][0]+points[0][0])/2, (points[len(points) - 1][1]+points[0][1])/2)   # 基准点选为第一个点和最后一个点连线边上的中点

    for i in range(0,len(points)-1):
        p1 = points[i]
        p2 = points[i + 1]

        db1 = coors_distance(pb[1], pb[0], p1[1], p1[0])   # 根据维度转化成经纬度距离
        d12 = coors_distance(p1[1], p1[0], p2[1], p2[0])
        d2b = coors_distance(p2[1], p2[0], pb[1], pb[0])

        hc = (db1+d12+d2b)/2   # db1是基准点和p1的距离，d12是p1和p2的距离，d2b是p2和基准点距离
        triArea = math.sqrt(hc * (hc - db1) * (hc - d12) * (hc - d2b))
        area += triArea

    return area
